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Abstract 


Using first-principles calculations, we have investigated the evolution of band-edges in 
few-layer phosphorene as a function of the number of P layers. Our results predict that mono- 
layer phosphorene is an indirect band gap semiconductor and its valence band edge is ex¬ 
tremely sensitive to strain. Its band gap could undergo an indirect-to-direct transition under a 
lattice expansion as small as 1 % along zigzag direction. A semi-empirical interlayer coupling 
model is proposed, which can well reproduce the evolution of valence band-edges obtained by 
first-principles calculations. We conclude that the interlayer coupling plays a dominated role in 
the evolution of the band-edges via decreasing both band gap and carrier effective masses with 
the increase of phosphorene thickness. A scrutiny of the orbital-decomposed band structure 
provides a better understanding of the upward shift of valence band maximum surpassing that 
of conduction band minimum. 


Recently, a novel two-dimensional material, few-layer black phosphorus (phosphorene) has been 
synthesized by micro-mechanical exfoliation technique. It arouses great interest of researchers 
due to many unparalleled properties superior to or not found in other members of the 2D materials 
family. For example, it exhibits a carrier mobility up to 1000 cm^/V-s and an on/off ratio up 
to 10"^ ~ 10^ for the phosphorene transistors at room temperature. Its peculiar puckered hon¬ 
eycomb structure leads to significant anisotropic electronic and optical properties on zigzag and 
armchair directions. Remarkably, the band gap of phosphorene is thickness-dependent, varying 
from 0.3 eV in the bulk limit to ~2.2 eV in a monolayer with a direct band gap character. 02,7,9-18 
Although the structural, electronic and optical properties of few-layer phosphorene have been 
extensively studied both experimentally and theoretically, there are issues that remain 

controversial or even unexplored, to the best of our knowledge. For example, there is a debate 
on whether single-layer phosphorene could have a direct-band gap.^^“^^ Cai et al. attributed the 
thickness-dependence of band gap to the quantum confinement effect. Nevertheless, they noted that 
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monolayer phosphorene is an exception. Tran et al. predicted that the decay of the band gap of 
phosphorene with increasing thickness is significantly slower than the usual quantum confinement 
result. This suggests that the quantum confinement effect alone is not enough to describe well 
the evolution of band gap as a function of the number of P layer. Another unresolved puzzle 
is that the upward shift of the valence band maximum (VBM) is faster than that of the conduction 
band minimum (CBM).^^’^^ Obviously, our understanding of the thickness effect on the evolution 
of the band-edges in phosphorene is still incomplete. 

In this Letter, we attempt to answer these questions by examining the fine structure of band 
edges in few-layer phosphorene using first-principles calculations based on the density functional 
theory (DFT).^*’^^ Our results demonstrate that monolayer phosphorene has an indirect band-gap 
and its valence band edge (VBE) is rather sensitive to the strain. A strain as small as 1% along 
zigzag direction is enough to induce an indirect-to-direct band gap transition in monolayer phos¬ 
phorene. Based on our first-principles results, we propose a semi-empirical interlayer coupling 
(SEIC) model to interpret the evolution of the valence band-edges as a function of the P number 
of layers in few-layer phosphorene. This simple model can reproduce fairly well the evolution of 
valence band-edges determined by first-principles calculations, a strong indication that interlayer 
interaction is primarily responsible for the evolution of the band-edges in few-layer phosphorene. 

Our total energy and electronic structure calculations were performed using the Vienna ab initio 
simulation package (VASP).^®’^^ The electron-ion interaction was described using projector aug¬ 
mented wave (PAW) method and the exchange and correlation were treated with generalized 
gradient approximation (GGA) in the Perdew Burke Ernzerhof (PBE) form.^"^ It is well known 
that the GGA method underestimates the semiconductor band gaps. However, in this study, we 
focus on the evolution of band-edges in few-layer phosphorene; the band dispersions calculated 
by GGA method exhibit similar features with the exception of the relative position of the valence 
and conduction bands when comparing with the other higher-level methods, such as hybrid DET 
or GW. We used a cutoff energy of 400 eV for the plane wave basis set, which yields total energies 
convergence better than 1 meV/atom. Previous DPT calculations have shown that the interlayer 
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van der Waals (vdW) interaction need to be considered for a proper description of the geometrical 
properties of black phosphorus. We therefore incorporated the vdW interactions by employ¬ 
ing a semi-empirical correction scheme of Grimme’s DFT-D2 method in the following calculations 
unless otherwise stated, which has been successful in describing the geometries of various layered 
materials. 

For the single-layer phosphorene, we also carried out local density approximation (LDA) us¬ 
ing the Ceperley-Alder exchange correlation potential as parametrized by Perdew and Zunger, 
quasiparticle GOWO and hybrid DFT^^“^^ calculations for the purpose of comparison. To achieve 
good convergence of dielectric function in the GOWO calculations, we used a large number of 320 
bands for the unit cell of monolayer. The converged eigenvalues and wavefunctions, as well as 
the equilibrium geometry obtained from PBE functional, were chosen as the initial input for the 
GOWO calculations. Note that in the GOWO calculations only the quasiparticle energies were re¬ 
calculated self-consistently in one iteration; the wave-functions were not updated but remain fixed 
at the PBE level. For visualization purpose, the GOWO band structure were interpolated to a finer 
grid using an interpolation based on Wannier orbitals as implemented in the WANNIER90 code.'^^ 
In the HSE06 approach, we here employed a revised scheme as proposed by Heyd, Scuseria, and 
Ernzerhof (HSEOb)."^^’'^'^ Additional computational settings and a more detailed discussion of the 
ground-state properties were given in our previous DET study of the native defect properties in 
phosphorene.^^ 

In the slab model of few-layer phosphorene, periodic slabs were separated by a vacuum layer of 
15 A in the c direction to avoid mirror interactions. A 8x6x1 k-mesh including F-point, generated 
according to the Monkhorst-Pack scheme,was applied to the Brillouin-zone integrations. On 
geometry optimization, both the shapes and internal structural parameters of pristine unit-cells 
were fully relaxed until the residual force on each atom is less than 0.01 eV/A. The fine structures 
of band structures were calculated by sampling 301 k-point along each high symmetry line in 
reciprocal space. The band structures and band edges were aligned with respect to the vacuum 
level which was determined by aligning the planar-averaged electrostatic potential in the vacuum 
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region far from phosphorene layer. This provides a deeper intuitional insight into the physieal 
origins of band edge evaluation. We used the VASPKIT eode for post-proeessing of the VASP 
ealeulated data."^^ 

We begin by investigating the global band strueture of monolayer phosphorene, aiming to find 
the loeations of both VBM and CBM in the two-dimensional reeiproeal lattiee spaee. In Fig. 1, we 
find that the band dispersion of VBE along the waveveetor ky (eorresponding to armchair direetion) 
near the high symmetry F point is rather signifieant, implying a small effeetive mass of the hole- 
earrier. By eontrast, the band dispersion of VBE along the waveveetor kx {zigzag direetion) is very 
flat, a signature of very heavy hole. A similar feature is also observed for the eonduetion band 
edge (CBE) dispersion. This means that highly anisotropie effeetive mass for both eleetron and 
hole earriers would be observed in monolayer phosphorene. ^ Meanwhile, it appears that both VBM 
and CBM are loeated at the E point in the Brillouin zone, seemingly indieating that the monolayer 
has a direet band gap. 



1 . 0 “' 


Eigure 1: (Color online) (a) Eayered strueture, (b) two-dimensional Brillouin zone, and (e) the 
highest valenee and lowest eonduetion bands of monolayer phosphorene. Green and red eurves are 
the projeetions of global band edges of monolayer onto the xz and yz planes. The energy level of 
vaeuum is set to zero. 

To gain further insights into the underlying physies, the orbital-projeeted band strueture of P 
atom (“fat bands”) are shown in Pig. 2. The band gap for monolayer phosphorene is predieted 
to be 0.91 eV by using PBE approaeh. The lowest-lying valenee bands are essentially 5 states. 
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characterized by large dispersion along F-X direction. Since the variation of kx has marginal effect 
on the py orbital, the bands consisting mainly of py states are relatively flat along F-X direction. 
The opposite is true for the bands with px character in F-Y direction. It is worth mentioning here 
that the Px and Py orbitals are not orthogonal to each other after hybridization. The VBE are mainly 
derived from the o-o bonding states between p^ orbitals (weight percentage: 90.7%) in different P 
sublayer. For the CBE, the contributions of s (14.8%) and py (19.6%) become comparable with p^ 
(65.6%). Clearly, the interlayer interaction has a stronger effect on the p^ orbitals than on others. 
It is therefore expected that the position of VBE is more sensitive to the number of P layers than 
that of CBE due to the interaction of p^ between neighboring P layers. 

A closer look of CBE and VBE are displayed in panels (b) and (c) respectively. It is seen 
clearly that for monolayer phosphorene the CBM is exactly located at F point, but the VBM shifts 
slightly away from F point by around ±6.33x 10“^ along the F-X direction. This means that 
monolayer phosphorene is not a direct-band gap semiconductor. This special position is labeled as 
A point and its energy and position deviation from F point are denoted as AEa and respectively. 
Based on the k ■ p perturbation theory, Ei et al. pointed that the indirect band gap character 
of monolayer phosphorene originates from the coupling among s, py and p^ orbitals since they 
have the same symmetry of Fj Our calculations show that the inclusion of spin-orbital coupling 
(SOC) shifts the VBE (CBE) toward (away from) Fermi level by around 8 meV (7 meV). However, 
the relevant physics discussed above is not affected qualitatively. Thus, the SOC effect will not be 
considered in the following discussion. 

For the purpose of comparison, we summarize in Table 1 the values of AEa and lSk\ calcu¬ 
lated using various methods with or without DFT-D2 correction. All methods predict that the 
monolayer phosphorene is not an exact direct gap semiconductor, in accordance with a recent GW 
work. Most noticeably of all, the magnitude of AE^ yielded by HSE06-D2 method even reaches 
~100 meV. We remind the reader that the calculated values of AEa listed in Table 1 correspond 
to different equilibrium lattice constants. Further HSE06-D2 calculations using the optimal lat¬ 
tice constants obtained with PBE-D2 (EDA) give AEa= 60 (93) meV and AkA=±0.15 (0.16) 
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Figure 2: (Color online) Orbital-projected band structure of monolayer phosphorene (a), and fine 
structures of CBE (b) and VBE (c) around the E point. The width of line in (a) indicates the weight 
of component. Blue solid and red dashed lines in (b) and (c) represent the band edges with and 
without spin-orbital coupling respectively. The vacuum level is set to zero. 


respectively. Clearly, the magnitude of AE^ is sensitive to the P-P bond length. Since the PBE- 
D2 method still overestimates the lattice constants for bulk black phosphorus,we speculate that 
the overestimation may extend to few-layer phosphorene. This might be a non-negligible factor 
to the very tiny value of AE/^. Moreover, EDA or GGA severely underestimates the semicon¬ 
ductor band gaps, the HSE06 method is expected to yield more accurate AEa for given lattice 
constants. In other words, the PBE method underestimates the energy difference between A and 
r when compared with HSE06 result. Interestingly, we find that with a vacuum thickness of 25 
A, GOWO method gives the values of AEa and Ak\ very close to those of EDA. It is, however, 
worth mentioning that the GOWO calculated band gap are dependent on vacuum thickness as we 
shall demonstrate below. Unfortunately, no experimental data are available for comparison with 
the predicted values at present. 

In view of the fact that GW calculations performed by different groups yield inconsistent band 
character of monolayer phophorene, we also performed GOWO calculations to examine the 
band structure of monolayer as a function of vacuum thickness. We find that the band struc¬ 
ture character of monolayer phosphorene, including both band dispersion and band gap, is rather 
sensitive to the vacuum thickness. As seen in Eig. 3 (a), the band dispersions of monolayer con- 
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Table 1: The calculated lattice constants a (A, zigzag direction) and b (A, armchair direction), 
P-P bond length d (A), energy and location deviation AE^ (meV) and AA^a (A ^) of A point 
away from F point in monolayer phosphorene, using various methods combined with/without 
DFT-D2 method. 


XC method a b d AFa AAa 


EDA 

3.27 

4.38 

2.20- 

'2.23 

12 

±0.13 

PBE 

3.30 

4.62 

2.22- 

'2.26 

0.6 

±0.06 

PBE-D2 

3.30 

4.57 

2.22- 

'2.25 

0.7 

±0.06 

HSE06-D2 

3.27 

4.53 

2.19- 

'2.22 

99 

±0.19 

GOWO 

3.30 

4.57 

2.22- 

'2.25 

11 

±0.12 
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Figure 3: (Color online) Band strueture of monolayer phosphorene given by GOWO (a) and fine 
struetures of CBE (b) and VBE (c) around the F point. Blue and red lines represent the results of 
a slab model with a vaeuum thickness of 25 A and 15 A respectively. The vacuum level is set to 
zero. 
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figurations with different vaeuum thiekness (25 A versus 15 A) are qualitatively similar. However, 
they have different band eharaeters when we look into the fine strueture of VBE in Fig. 3(o). More 
speeifieally, with a vaeuum thiekness of 25 A, GOWO prediets an indireet band gap, in agreement 
with our PBE and HSE06 results. If we adopt a vaeuum with a thiekness of 15 A in the slab model, 
GOWO prediets a direet band gap. This sensitivity might be the reason for the diserepaney in previ¬ 
ous GOWO results reported in literature. In addition, the gap value also varies linearly as a funetion 
of the inverse of vaeuum thiekness (see Fig. 4). Similar behavior was also found in other two 
dimensional materials,refieeting the long-range nature of the self-energy 'L=iGW.^^ Finally, 
the asymptotie gap value of 2.01 eV for monolayer phosphorene in the limiting ease of infinite 
vaeuum is obtained by an extrapolation seheme, very elose to the reeently reported experimental 
value of 2.2 eV.^ 



0.02 0.04 0.06 0.08 

Inverse of vacuum thickness (A^ ) 


Figure 4: (Color online) Band strueture of monolayer phosphorene as a funetion of the inverse of 
vaeuum thiekness given by GOWO. 


Next, we investigate the infiuenee of lattiee eonstants on the AFa by applying uniform out- 
of-plane (e^) or in-plane strain to the unit eell of mononlayer. The latter ineludes three eases, 
namely, the uniaxial strain along x (zigzag) or y (armehair), and the biaxial strain along both x and 
y direetions, denoted by e^, Ey and £j^y respeetively (see Fig. 1). Taking the ease of the uniaxial 
strain along zigzag direetion as an example, the applied strain is defined as where Ux and 

^A'O 

UxQ are the lattiee eonstants along the x direetion for the strained and relaxed struetures. A positive 
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(negative) value of e corresponds to a tensile (compressive) strain. 



Figure 5: (Color online) Evolution of VBE for monolayer phosphorene in response to strains 
determined by PBE-D2. For explanation of £x. By, £z and £xy, please see text. 

As seen in Fig. 5, a uniform out-of-plane or in-plane compression lowers the eigenvalue at 
the r point. This leads to more significant indirect-gap character of monolayer phosphorene. In 
contrast, the system undergoes a transition from indirect to direct band gap when tensile strains are 
applied. Our calculations predict that a small tensile strain of around 1% along zigzag direction 
will induce such transition. This implies that the band gap character of thin phosphorene films 
grown on different substrate could be different due to the introduction of strain by lattice mis¬ 
match. In contrast, the evolution of VBE is not sensitively dependent on the strain along armchair 
direction. It is readily understandable if we consider the fact that the VBM is located at the A point 
in F-X path (corresponding to zigzag direction) as pointed out above. Previous first-principles 
calculations also showed that both the electronic and optical properties of single-layer phospho¬ 
rene depend strongly on the applied strain. To gain insight into such observed trend in the 
VBE of monolayer phosphorene, we examine the bonding characters of F point by plotting the 
band-decomposed electron density in Fig. 6. 

We find that the depth of concave-like shape of VBE near F point is associated with the strength 
of ppo bonding states between p^ orbitals of P atoms. More specifically, stronger bonding strength 
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between orbitals induces more significant indirect band gap character. Bonding strength can 
be enhanced by reducing the P-P bond length in different sublayers, i.e., applying a compressive 
out-of-plane or in-plane strain. Phosphorus atoms have five electrons on 3p orbitals, three of 
which participate in the formation of three covalent a-bonds with three neighboring P atoms by 
sp^ hybridization in a puckered honeycomb structure, and the remaining two electrons occupying 
a lone pair orbital oriented out-of-plane. Note that the magnitude of is more sensitive to in¬ 
plane strain than out-of-plane one. This means that the indirect-gap character of VBE is also likely 
associated with the coulomb repulsion between electrons from the lone pair orbitals of neighboring 
P atoms. 



Figure 6: (Color online) Charge density plots of the F point when monolayer phosphorene is 
subject to varies strains. The isosurface is 0.01 e/Bohr^. 

When two or more isolated monolayers move close to each other to form multilayer phospho¬ 
rene, the corresponding degenerated energy levels of different monolayers become non-degenerated 
due to the interlayer coupling, as is shown in Fig. 7. This leads to the splitting of the correspond¬ 
ing bands, which in turn push VBM and CBM to higher and lower energies respectively. As a 
result, the magnitude of band gap decreases with increasing the number of layers. On the basis 
of the careful analysis of the fine structure of conduction and valence bands, as shown in panels 
(b) and (c), we conclude that bilayer phosphorene is a direct gap semiconductor. In other words, a 
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transition from indirect to direct band gap occurres when going from monolayer to bilayer. 



Figure 7: (Color online) Orbital-projected band structure of bilayer phosphorene (a) and fine struc¬ 
tures of CBE (b) and VBE (c) around the E point. The width of line in (a) indicates the weight 
of component. Blue solid and red dashed lines in (b) and (c) represent the band edges with and 
without spin-orbital coupling respectively. The vacuum level is set to zero. 


Eigures 8 and 9 show the evolution of band-edges along E-X and E-Y paths as a function 
of the number of P layers respectively. One can find the band dispersion of VBE along E-X 
path becomes more significant with increasing film thickness due to the strong effect of interlayer 
coupling on orbital, suggesting that a significant decrease of hole-carrier mass is to be observed 
with increasing the number of layers. Generally, a decrease in the carrier effective mass means 
an enhanced mobility. Therefore, few-layer phosphorene is likely to have better performance than 
monolayer phosphorene in electronic devices. 

We find that the band dispersions of both VBE and CBE along E-Y path change from parabolic¬ 
like shape for monolayer to (nearly) linear-like shape for quadrilayer, suggesting that the Dirac- 
point observed in graphene could be also realized in few-layer phosphorene by tuning its band 
gap value to zero eV. On the other hand, the VBE along the E-X direction remain parabolic-like 
shape, indicating a significantly anisotropic band structure of quadrilayer. We note a very recent 
experimental study has reported that the tunable band gap and anisotropic Dirac semimetal state 
can be achieved by sprinkling potassium atoms on top of multilayer phosphorene.^^ Several recent 
theoretical studies also predicted that the Dirac cone in phosphorene can be achieved by application 
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Figure 8: (Color online) The evolution of the band-edges along F-X direction as a function of the 
number of layers calculated by DFT and SEIC model, (a) monolayer, (b) bilayer, (c) trilayer, (d) 
quadrilayer. The vacuum level is set to zero. 
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Figure 9: (Color online) The evolution of the band-edges along F-Y direction as a function of the 
number of layers calculated by DFT and SEIC model, (a) monolayer, (b) bilayer, (c) trilayer, (d) 
quadrilayer. The vacuum level is set to zero. 
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of an external eleetrie field or strain. 

For quadrilayer, it should be pointed that the PBE method prediets a tiny gap value of around 
0.1 eV, whieh is far smaller than the value of 0.71 (1.08) eV given by HSE06 (GWO).^^ Therefore, 
we attribute the linear dispersion of quadrilayer to the strong interaetion between valenee band and 
eonduetion band as a result of the severe gap underestimation by PBE. Further HSE06 ealeulations 
reveal that the band dispersion of quadrilayer remains a parabolie-like shape but has a tendeney 
toward linear-dispersion when going from monolayer to the bulk limit, as shown in Supporting 
Information Figure SI. It is found that the VBM shifts toward Fermi level by 0.7 eV while the 
CBM shifts downward by 0.1 eV when going from monolayer to quadrilayer. The faet that the 
weight of Pz in VBE is larger than that in CBE is likely to be responsible for sueh observations. A 
similar trend has been found for M 0 S 2 . The global evolution of both VBE and CBE for bilayer, 
trilayer and quadrilayer systems is depieted in Supporting Information Figure S2. 

It is now elear that the behavior of VBE is erueial for the eleetronie properties of phosphorene 
systems. Thus, we take the ease of VBE as an example and propose a SEIC model to deseribe 
the evolution of VBE as a funetion of the number of P layers. As an extended Hiiekel moleeular 
orbital theory,the SEIC model ean deseribe the interlayer interaetion in layered materials. In 
this model, a multi-layer phosphorene is seen as a “moleeule”, a monolayer is taken as an “atom”, 
the energy bands in a monolayer are analogue of eleetron orbitals, and the interlayer eoupling is 
the hybridization of the orbitals originating from different “atom” with the same symmetry.. Based 
on k ■ p perturbation theory, the band dispersion of VBE near the F point in monolayer ean be 
approximately expressed as 


E\vikxi ky') — ky') — £0 T (^xk^ T Pxk^ 

+ (Xyky + Pyky, 


( 1 ) 


where £o=-4.92 eV is the eigenvalue of F point of monolayer. The parameters ttx, Oiy, jix and are 
0.04 eV-A^, -30.0 eV-A^, -80.0 eV-A"^ and -100.0 eV-A^ respeetively, indieating the anisotropie 
eleetronie properties on zigzag and armehair direetions. These values ean be obtained by fitting 
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the DFT calculated band dispersion of VBE for monolayer phosphorene. Here kx and ky are given 
in units of Inla and Inib respectively, where a and b are the corresponding equilibrium lattice 
constants of few-layer phosphorene systems. Since we focus only on the interpretation for the 
evolution of the valence band-edges along zigzag {ky= 0 ) and armchair (kx= 0 ) directions of few- 
layer phosphorene, the contribution of coupling between kx and ky to the valence band dispersions 
is neglected. 

As for the multilayer, the interlayer interaction is described by a coupling term of 7v(k). For 
simplicity, the interaction between non-neighboring layers is neglected without loss of accuracy. 
In the case of VBE, the coupling of p^-Pz orbitals in nearest-neighboring layers is only considered 
since the VBE is mainly derived from pz state as pointed above. Then the interlayer coupling term 
of VBE is defined as 

Jvij^xiky) = flVxkx T ^ykyf (2) 


where /i, Vx and Vy are the coupling parameters. 

The Hamiltonian of valence band dispersions near the F point in bilayer is 


H^lvikxi ky) — 


H\v Jv 

Jv ki\v 

V / 


and the eigenvalues of Eqs. (??) are 


( 3 ) 


^Iv^k-X-t ky) - E\x(kxt ky) i7v(^X; ky) . 


( 4 ) 


Based on the DET-calculated results of valence band dispersions for bilayer as displayed in 
Eig. 8 (b), the parameters p, Vx and Vy are 0.44 eV, -4.18 eV-A? and 6.0 eV-A^ respectively. 
One can find that the energy of the first and second highest valence band of bilayer are shifted 
upward and downward respectively by Jy(kx,ky) eV, with respect to monolayer. Note that Jx(kx,ky) 
is wavevector-dependent and reaches its maximum /i=0.44 eV at the F point. Clearly, our SEIC 
model reproduces the indirect-direct band gap transition when going from monolayer to bilayer 
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as predicted by our DFT results. Next, we take trilayer as an example to study the role of inter¬ 
layer coupling on the evolution of VBE in multilayer system. The Hamiltonian of trilayer and the 
corresponding eigenvalues can be described as 


and 


//iv Jv 0 


ky) — 


Jv H 


Iv 


Jv 


0 Jv Hiv 


— Eivikxiky) \/2Jv{kx,ky^, 
^’iv^^xiky) =Eiv{kx,ky^, 

El^ikx.ky) = Eivikx.ky) -f V2Jv{kx,ky). 


( 5 ) 


( 6 ) 


The analysis is similar for the systems with more layers. For the n-layer system, the evolution 
of its VBE is determined by Env{kx, ky)=Eiv+'2Jv-cos-^n, where k=l, 2,..., n is the index number 
of splitting valence band. In the bulk limit (n —oo), the eigenvalue of E point is expected to be 
Eiv+27v=- 3.16 eV with respect to the vacuum level. In other words, compared with monolayer the 
gap value of bulk black phosphorus reduces at least 0.88 eV if one also consider the evolution of 
CBE toward Fermi level. This prediction agrees well with our previous DFT-PBE results. 

Overall, the SEIC model results are in good agreement with those obtained from DFT, as shown 
in Figs. 8 and 9. However, the SEIC model predicts that the VBE of quadrilayer or more-layer 
systems along E-Y direction still remain quadratic-like shape, instead of showing linear-like shape 
from DFT calculations [see Figs. 9 (c) and (d)]. This disagreement is probably caused by the 
negligence of the stronger coupling between VBE and CBE when they get close in the quadrilayer 
or thicker systems. However, it should be emphasized here that the artificial linear dispersion 
of quadrilayer is derived from the problem of the band gap underestimation related to the local 
or semilocal approximations of the exchange correlation functional in DFT. Based on the first- 
order k ■ p perturbation theory, we use a degenerate two-band model to describe the coupling 
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between the valenee and conduetion bands for these systems. Our results indicate that this coupling 
does lead to the linear-dominated dispersion along armchair direction and the quadratic-dominated 
dispersion along zigzag direction respectively. 

We expect our SEIC model is also applicable to the evolution of CBE. The evolution of CBE is 
different from that of VBE. As an illustrative example we choose the case of trilayer system, one 
can observe in Eig. 8 (c) that the energies of the second and third conduction bands near the E point 
become very close to each other; while the second and third valence bands is separated by an energy 
of y/2Jy{kx,ky) eV. This means the coupling among 5, py and become noticeable and cannot 
be neglected. Eurther calculations show that the relative positions of the first four conduction 
bands at the F point can be qualitatively determined if we include the s and py states into the 
SEIC model to describe the evolution of the CBE. Previous theoretical studies suggested that the 
decreasing tendency of band gap with increasing film thickness is likely an outcome of the quantum 
confinement effect. We argue that it would be sufficient to explain such an observed tendency 
using a simple SEIC model. Clearly, the dramatic decrease of hole mass is also attributed to 
interlayer coupling. We hope that this systematic study serves as a guideline for the interpretation 
of the evolution of band edges in other two-dimensional layered materials. 

To conclude, through first-principles calculations and a semi-empirical interlayer coupling 
model, we have uncovered the nature and evolution of the valence band-edge as a function of 
the number of P layers in few-layer phosphorene. Our results predict that the interlayer coupling 
play a vital role in determining the decreasing trend for both band gap and carrier effective masses 
with the increase of phosphorene thickness. Also, the interlayer coupling leads to a transition from 
indirect to direct when changing from monolayer to few-layer phosphorene. The analysis of the 
orbital-decomposed band structure reveal that the upward shift of valence band maximum is faster 
than that of conduction band minimum, due to the larger contribution of p^ in the former. 
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